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We explore the link between soft vibrational modes and local relaxation events in polymer glasses 
during physical aging, active deformation at constant strain rate, and subsequent recovery. A 
softness field is constructed out of the superposition of the amplitudes of the lowest energy normal 
modes, and found to predict up to 70% of the rearrangements. Overlap between softness and 
rearrangements increases logarithmically during aging and recovery phases as energy barriers rise 
due to physical aging, while yielding rapidly rejuvenates the overlap to that of a freshly prepared 
glass. In the strain hardening regime, correlations rise for uniaxial tensile deformation but not for 
simple shear. These trends can be explained by considering the differing degrees of localization of 
the soft modes in the two deformation protocols. 

PACS numbers: 64.70.P, 81.05.Kf, 61.43.Fs, 61.43.Bn 


I. INTRODUCTION 

The mechanical behavior of solids is often successfully 
described by continuum mechanics of a homogeneous 
medium. Polymer glasses, however, exhibit structural 
heterogeneity at the nanoscale. For instance, local elastic 
shear moduli evaluated in small volume elements contain¬ 
ing several tens of monomers can differ substantially from 
the bulk macroscopic value [1]. As a result, amorphous 
materials can be separated into stiff and soft regions. At 
the same time, glassy dynamics at the monomer level is 
heterogeneous, with groups of particles relaxing rapidly 
in a collective manner while other regions remain essen¬ 
tially immobile mm. Understanding the existence or 
absence of a link between local structure and dynamics 
is one of the important open questions in the physics of 
disordered materials. 

Recent molecular simulations of a variety of model 
glass formers have provided strong evidence that the lo¬ 
cations of structural rearrangements correlate with mea¬ 
sures of structural heterogeneity. For instance, simula- 
tons of a model polymer glass under tensile deformation 
show that nonafiine residual monomer displacements are 
largest in regions with small positive shear moduli [4]. 
Simulation studies of the athermal deformation of amor¬ 
phous two-dimensional binary mixtures also indicate that 
the nonafiine displacement field resulting from the macro¬ 
scopic shear deformations of a few percent is directly re¬ 
lated to the spatial structure of the elastic moduli mis], 
while little correlation is found to other structural vari¬ 
ables such as local density. 

An alternative description of structural heterogeneity 
that also untilizes a harmonic description has been devel¬ 
oped around the notion of soft vibrational modes. Amor¬ 
phous solids generally exhibit an excess amount of low 
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energy modes (Boson peak), and a large number of them 
are quasilocalized, i.e. they involve only very few par¬ 
ticles [ 7 ]. The particles involved in these modes exhibit 
structural differences in their neighborhood configuration 
[8] , suggesting that they might play the role of defects in 
crystals. Recently, it has become possible to identify the 
structural signatures of rearranging particles with gener¬ 
alized structure functions that describe radial and bond 
orientations [9]. A “softness field” formed by the super¬ 
position of the amplitudes of the lowest frequency normal 
modes also yields a heterogeneous partitioning into hard 
and soft regions. Multiple recent simulation studies in 
supercooled liquids pd>HT2] and amorphous packings un¬ 
der shear at zero [IMS] or finite temperature [16] con¬ 
firm that regions of large vibrational amplitude overlap 
with the loci of structural rearrangements. An advan¬ 
tage of the soft mode description over the elastic moduli 
approach is that the soft modes also anticipate the direc¬ 
tion of particle motion in a rearrangement UlllIH]. This 
correlation is robust and insensitive to the specific model 
system or diagnostic of the rearrangement. Low energy 
sound waves therefore scatter off flow defects and provide 
information about those collective particle motions that 
are most easily excited. 

In a recent contribution, we explored quantitatively 
the robustness of this softness field for the prediction of 
monomer relaxation in aging polymer glasses [17] for dif¬ 
ferent temperatures and aging times. We showed that 
rearrangements identified as rapid changes in the parti¬ 
cle positions (hops out of local cages) are up to seven 
times more likely at the softest regions, where the direc¬ 
tion of motion also aligns near perfectly with the local 
polarization vectors. These correlations are only erased 
once more than 50% of the entire system has undergone 
rearrangements. 

In the present work, we report a study of the correla¬ 
tion between irreversible rearrangements and local soft¬ 
ness in actively deformed polymer glasses. This is im¬ 
portant because previous studies have focused only on 
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steadily sheared two-dimensional monomeric glasses for¬ 
mers m m HI], while polymers are distinctly three- 
dimensional and exhibit strain hardening. Here we sub¬ 
ject the polymer to uniaxial tensile deformation (exten- 
sional flow) and compare to the simple shear protocol 
used in previous studies. Moreover, we explore subse¬ 
quently correlations between soft regions and local re¬ 
laxation events in the regimes of aging, deformation, 
and recovery. During the aging phase, correlations grow 
due to a descent into deeper energy minima, but return 
quickly to their as-quenched values once plastic flow sets 
in. In the uniaxial tensile deformation protocol, correla¬ 
tions rise again with increasing strain, which we show to 
arise from an increase in the degree of localization of the 
soft modes as the aspect ratio of the sample changes. In 
simple shear by contrast, this effect is absent and correla¬ 
tions remain near constant during deformation. During 
the recovery phase, correlations increase only weakly as 
inherent structure energies return to their undeformed 
values. 


II. METHODS 

We study the polymer glass using molecular dynamics 
techniques and the well-known finitely extensible non¬ 
linear elastic (FENE) bead-spring model [19], which is 
an excellent computer glass-former Eoiini- Each linear 
polymer consists of 50 identical monomers with covalent 
bonds along the polymer backbone that were modeled 
by a non-linear stiff spring-like interaction. Inter- and 
intra-chain interactions between not covalently bonded 
monomers are modeled with a 6-12 Lennard-Jones (LJ) 
potential, which was truncated and force-shifted [22] with 
a cut-off of 2.5 times the bead diameter for computa¬ 
tional efficiency, and to ensure that the Hessian is de¬ 
fined. The latter is important for the analysis of the 
vibrational spectrum. All quantities are stated in the 
usual LJ-units based on well energy particle d iameter 
a, mass m and characteristic time scale tlj = yjmo? ju^ 
which is just below twice the mean collision time. The 
simulation time step was set to At = 0.0075rLj, and all 
results are averaged over 20 independent simulation runs. 

We simulate N=10,000 particles comprising 200 poly¬ 
mers in an initially cubic simulation box with periodic 
boundary conditions. Our simulation protocol is a direct 
extension of our earlier study that quantified the corre¬ 
lation between soft modes and irreversible particle rear¬ 
rangements called hops in quiescent polymer glasses HZ]. 
The glass is generated by rapidly quenching an equili¬ 
brated melt with density p = 1.043 at constant volume 
(NVT) from T = 1.2 to T = 0.3 with a constant quench 
rate of T = 6.7 • 10“^. Eor the remainder of the simula¬ 
tion, the temperature is held constant. In a first step, the 
glass is aged without deformation at zero pressure (NPT) 
to the age tage = 7, 500, 000 at T = 0.3, which is just be¬ 
low the glass transition temperature. In the uniaxial ten¬ 
sile deformation protocol, the glass is then deformed in 


one direction with a constant strain rate e = 10“^, while 
the pressure perpendicular to the deformation axis is kept 
at zero. As a result, the simulation box shape changes 
from cubic to rectangular. The deformation ends at a 
final engineering strain e(t) = [L^it) — Lz{0)]/Lz{0) = 4 
with Lz being the simulation box size along the defor¬ 
mation axis. In a final step, both deformation and baro- 
stat were turned off so that the system recovers at fixed 
volume and the stresses relax. This post-deformation 
recovery regime was explored for a time tr = 150, 000. 
The simple shear protocol imposes a constant shear rate 
of the same magnitude by deforming the simulation box 
into a parallelepiped. Here no additional barostating is 
performed. 

Two key measurements are performed during the sim¬ 
ulation: At different times during aging, deformation 
and recovery, we first store a configuration snapshot for 
the identification of the soft modes and then detect irre¬ 
versible particle rearrangements in the whole simulation 
box immediately following the time of the snapshot. As 
introduced in full detail in Ref. HU, we calculate a soft¬ 
ness measure for each particle based on its participation 
in the low energy vibrational modes. This calculation 
was done in multiple steps: Starting from a snapshot, we 
first found the inherent structure using a combination of 
gradient descent and damped dynamics (EIRE [23]) al¬ 
gorithms with a minimal total force criterion. Then, the 
Hessian 
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is calculated from the inherent structure particle lo¬ 
cations r as well as the potential energy I/, and the 
Nm = 600 lowest energy eigenmodes are calculated us¬ 
ing ARPACK. The softness of a particle is defined as the 
superposition of the participation fractions in the low en¬ 
ergy vibrational modes Hnuniis], 
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Here, the polarization vector is the projection of the 
eigenvector of mode j on the degrees of freedom of parti¬ 
cle i. The softness field (j) depends on a single parameter, 
the number of low energy modes and the scaling 

factor is added to make the softness an intensive quan¬ 
tity in terms of An optimal value of = 600 (2% 
of the modes) is chosen based on our previous analysis 
HZ]. A particle i is ’’softer” the larger 0^, which is used 
to rank the particles according to their relative softness. 
The absolute value of (j) is not in itself meaningful, since 
the participation fractions are normalized quantities, i.e., 

Ef=i|efP = l. 

Eigure [^a-c) visualizes three exemplary low energy 
eigenmodes calculated at the end of the aging period by 
plotting their polarization vector fields. One can see that 
the modes (a) and (b) are quasi-localized, since large po¬ 
larization vectors are concentrated on a small number of 
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Figure 1. (a-c) Three exemplary lowest energy vibrational 

modes, visualized by their polarization vector fields (coloring 
indicates depth). The participation ratios as defined in Eq. @ 
are 0.12(a), 0.29(b) and 0.62(c). (d) Exemplary simulation 
snapshot with particle softness indicated by color (blue to 
red - small to large 0). The right side shows only the 10% 
softest particles. 


particles that are spatially clustered. The spectrum also 
exhibits extended, planar-wave like modes and an exam¬ 
ple is shown in panel (c). Panel (d) shows a snapshot of 
the system where each particle is colored according to its 
softness (blue is small (j) and red is large (j)). On the right 
side of the simulation box, we only show the 10% softest 
particles, illustrating the heterogeneous spatial distribu¬ 
tion. 

Irreversible particle rearrangements are identified us¬ 
ing a detection algorithm that was introduced in a pre¬ 
vious study [24], where it is explained in greater detail. 
In short, we monitor the trajectory of each particle on- 
the-fly during the simulation and isolate rapid changes 
or hops of particles by calculating the mean distance 
squared 



between the earlier (A) and later (B) half of the trajec¬ 
tory of particle i in a time window that moves with the 
simulation. Here, the averages {.)a [(•)b] are taken over 
all trajectory points in A [B] and [f^] is the mean 
position in the respective trajectory segment. A hop is 
detected when Phop > Pth = 0.21, and we record particle 
identifier, time of hop, as well as initial and final position. 
The value of the threshold is related to the plateau in the 
mean square displacement [24] . To isolate the irreversible 
rearrangements we exclude back-and-forth hops of a par¬ 
ticle between the same two positions. This is done by 
removing a sequence of two hops of the same particle, if 
the final position of the second hop is within a distance 
of y/Pth /2 of the initial position of the first hop. 


III. RESULTS 

We discuss the spatial correlation between softness 
field and monomer hops in terms of a predictive success 
rate 0 m- Here, the softness field is binarized into a soft 
spot map denoted by assigning a softness of one to 
the fraction / of particles with largest softness and zero 
to the other particles. We then calculate the fraction of 
the first = 100 hopping particles that are part of a 
soft spot, or 
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with hi = 1 if particle i is one of the first particles 
to hop after the measurement of the softness field, and 
hi = 0 otherwise. 

Figure]^ a) shows the predictive success rate measured 
at three different strains during the uniaxial tensile de¬ 
formation (qualitatively similar curves are found under 
simple shear). Stress-strain curves for the two deforma¬ 
tion modes are shown in Fig. [^b) and reveal post-yield 
strain hardening in both cases, albeit to a much lower 
degree in the simple shear case. It can be seen that the 
three strains shown in panel (a) correspond to the onset 
of deformation at the end of the aging period, right af¬ 
ter the yield peak, and the end of the strain hardening 
regime. The dotted line indicates 0 for randomly dis¬ 
tributed soft spots (no correlation) for comparison. For 
all strains, we find a larger than random correlation be¬ 
tween softness field and the occurrence of hops. 

To facilitate the analysis of trends in the spatial cor¬ 
relation during aging, deformation and recovery, we now 
focus on the predictive success rate at a coverage fraction 
of 30%. This is the fraction where the difference between 
measured 0 and uncorrelated value (dashed line) is max¬ 
imal. In Fig.j^d) we show the evolution of 0(30%) dur¬ 
ing the initial aging period. In agreement with the results 
reported in Ref. m, the correlation increases from 0.61 
to 0.71 as the age grows by four orders of magnitude. 
Panel (e) shows the predictive success rate measured at 
different points during uniaxial tensile and simple shear 
deformation. In the tensile case (solid line), the spatial 
correlation in the elastic regime at e = 0 is roughly equal 
to the value found in the quiescent state immediately 
prior to the deformation. At the yield strain e = 0.04, 
the correlation has decreased to 0.64 and at the end of the 
strain softening regime (e = 0.1) 0(30%) has reached the 
pre-aging value 0.61. Even larger decreases in fact below 
the pre-aging value are seen in the pure shear case. This 
reversal of the aging effects is consistent with mechani¬ 
cal rejuvenation, by which deformation strains exceeding 
the yield strain erase the thermal history and return the 
material to a freshly quenched glassy state m- 

Interestingly, trends in spatial correlation begin to di¬ 
verge between the two deformation modes in the strain 
hardening regime. Under uniaxial tension, we find that 
0(30%) monotonically increases with growing strain. 
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Figure 2. (a) Fraction of hops in soft spots as function of coverage fraction of soft spots measured in three deformation regimes 
for uniaxial tension: elastic (e = 0.0), strain softening (e = 0.1), and strain hardening (e = 4.0). The dotted line indicates 
no correlation, (b) Stress-strain curves for uniaxial tensile (solid line) and simple shear (dashed line). Vertical colored lines 
indicate the investigated deformations at strains e = 0.0, 0.04, 0.1, 0.5, 1.0, 2.0, 3.0, 4.0 and the inset shows the region of the 
yield point in more detail, (c) Stress recovery from uniaxial tensile deformation. The evolution of the predictive success rate 
reached at / = 0.3 is shown during aging (d), deformation (e) and in the recovery regime (f). 


The increase accelerates at large strains e > 2, reach¬ 
ing a value of 0.72 at e = 4, which is above the predictive 
success rate measured in the quiescent state prior to the 
deformation. By contrast, the correlation during simple 
shear saturates at the value of the freshly quenched glass 
for strains larger than one. In the recovery regime, in¬ 
vestigated only for the tensile case (panel (f)), 0(30%) 
first drops slightly but then increases again weakly as the 
material continues to age. 

What is the origin of the different behavior of the 
softness-rearrangement correlation in tensile vs. simple 
shear deformation? In order to answer this question, we 
analyze the low energy vibrational mode spectrum dur¬ 
ing the deformation period more carefully. The extent of 
localization of a mode j can be quantified by the partic¬ 
ipation ratio 


(Ef=i(ef)^)' 


( 4 ) 


A value of Pj = 1 means that all particles are participat¬ 
ing equally in mode j, whereas a small value indicates 
that the mode is quasi-localized around a few active par¬ 
ticles. In Fig. |^(a) we show the participation ratio P{uj) 


as function of mode frequency uj. Since the vibrational 
spectrum is a feature of the inherent structure, P{oj) at 
strain e = 0 reflects the quiescent state at the end of 
the aging period. The participation ratios at the end 
of the strain softening regime (e = 0.1) are nearly un¬ 
changed, with a slight shift of the extended modes (large 
participation ratio) towards smaller frequencies for both 
deformation modes. 

Important differences emerge in the strain hardening 
regime at e = 4. While P{uj) for simple shear is un¬ 
changed, uniaxial tension deformation results in two new 
features: First, we find modes with large participation 
ratio at much smaller frequencies. This change is due to 
the large (400%) elongation of the simulation box along 
the deformation axis in our uniaxial tension simulation, 
which allows modes with larger wavelength to ”fit” into 
the simulation volume. In general one expects extended 
modes to emerge on the low energy side of the spectrum 
as the system dimension increases, since the material 
looks increasingly homogeneous on large length scales. 
At same time, however, we find that the participation 
ratios of the modes at cj > 1.5 are reduced compared to 
the undeformed system. With increasing linear dimen¬ 
sion, the likelihood of finding a low energy quasilocalized 
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Figure 3. (a) Mean participation ratio as function of eigen- 
frequency at three values of strain for uniaxial tensile (o) and 
simple shear (v) deformations., (b) Mean participation ra¬ 
tion of all Nm modes used for the softness field calculation as 
function of total strain. 

mode is therefore enhanced. This effect is absent in the 
simple shear case, during which the simulation box peri¬ 
ods do not change. 

In fig. I^b) we show the evolution of the mean partic¬ 
ipation ratio (P(co’)) of all = 600 modes used for 
the calculation of the softness field with deformation. 
For tensile deformation, we find a gradual reduction of 
the average participation ratio, indicating that the soft 
modes become more and more localized as the deforma¬ 
tion grows. By contrast, the average participation ratio 
does not change during simple shear. These results sug¬ 
gest that the change in participation ratio is merely a 
consequence of the changing box dimensions during uni¬ 
axial strain. In order to separate deformation effects from 
geometry, we have equilibrated the same polymer glass in 
a simulation box that has identical shape to a deformed 
polymer at 400% strain and computed the distribution 
of participation ratio at zero load. These distributions 
are identical to those found at the end of the hardening 
regime, which confirms that the change in {P{uj)) is of 
geometrical origin. 

We now show that the increasing degree of localization 
of the soft modes under uniaxial deformation explains 
the rise of the overlap 0(30%) between softness and re¬ 
arrangements. To this end, we compare the evolution 
of 0(30%) to changes in the mean participation ratio in 


the top panels of Fig. To facilitate the comparison, 
we introduce a scaling that normalizes these quantities 
by the extremal change during the simulation run. The 
relative change of a quantity Q with regard to a refer¬ 
ence value Qo from the freshly quenched state is defined 
as {Q - Qo)/(Qext - Qo), where Qext is the extremal 
value measured during the aging, deformation and in the 
recovery regime. For the spatial correlation we use the 
maximal value, while for the mean participation ratio we 
use the minimal value. The top left panel shows that the 
mean participation ratio does not change during aging, 
while the overlap 0(30%) rises roughly logarithmically. 
The overlap between hops and soft modes improves as 
the system reaches deeper energy wells and the energy 
barriers for rearrangement rise with age. During uniax¬ 
ial tensile deformation (middle panel), we find that after 
the initial rejuvenation drop, 0(30%) rises in lockstep 
with the mean participation ratio (P(co’)). This obser¬ 
vation suggests that the improving overlap of 0(30%) 
in the hardening regime is indeed due to an increased 
amount of highly localised modes. During simple shear 
by contrast, neither overlap nor mean participation ra¬ 
tio increase with strain. (P(cj)) is again constant in the 
recovery regime where the simulation box size does not 
change. 

The bottom panels of Fig. [^provide as additional in¬ 
formation the evolution of the inherent structure energy, 
separated into contributions from pairwise LJ interac¬ 
tions and covalent bonds. During aging (left panel), we 
see the well-known logarithmic aging increase of 
while the mean bond energy remains constant. 

These trends mirror those in the panel above for the 
overlap. In the deformation stage (middle panel), 
decreases as the glass is pulled up higher on the energy 
landscape, while increases as more and more en¬ 

ergy is stored in the covalent bonds of the aligning poly¬ 
mer chains. While this increase of with strain is 

also broadly consistent with the rise of the overlap mea¬ 
sure in the uniaxial tensile deformation, it disagrees in 
simple shear deformation where the overlap is constant. 
Moreover, in the recovery regime (right panel), 
drops rapidly while 0(30%) is very weakly increasing. 
These observations provide further evidence that strain 
hardening does not influence the softness map - rear¬ 
rangement correlations. Our recovery regime is just long 
enough for both inherent structure energies to reach their 
undeformed values. 


IV. DISCUSSION 

We have investigated the correlation between local re¬ 
laxation events, identified as irreversible monomer hops 
out of local cages, and a softness map constructed from 
the superposition of low energy vibrational modes, in 
polymer glasses during physical aging, deformation and 
structural recovery. In the as-quenched glass, hops are 
about twice as likely to occur on soft spots than else- 
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Figure 4. Comparison of spatial correlation, participation ratio (top row) and potential energy stored in the LJ interaction and in 
the covalent bonds (bottom row) during aging, deformation and recovery. The relative change is defined as (Q —Qo)/ (Qext — Qo), 
where Q is a placeholder for one of the measured quantities (see legend), and Qext is an extremal value observed during the 
simulation run (maximum for 0(30%), and minimum for {P(uj)), U^s). 


where in the polymer. This correlation increases loga¬ 
rithmically to about 2.4 over 4 decades of aging time and 
rises because physical aging amplifies the influence of the 
potential energy landscape on the dynamics. Both simple 
shear and uniaxial tensile deformation rapidly rejuvenate 
this correlation once plastic flow sets in. Further uniax¬ 
ial tensile deformation increases the soft spot-relaxation 
overlap again due to an increase in quasilocalized modes 
as the sample elongates in the tensile direction, while 
no such effect occurs in simple shear. Strain hardening, 
although present in our model polymer, does not seem 
to play a major role in these trends. Accordingly, the 
present effects can be expected in short molecule glass 
formers as well. 

The present results indicate that two factors control 
the ability of soft modes to predict structural rearrange¬ 
ment: position on the energy landscape, and degree of 
localization of the lowest energy modes. The lower the 
inherent structure energies and the higher the barriers. 


the better the ability of the potential energy landscape (a 
T=0 property) to predict dynamics at elevated temper¬ 
ature near the glass transition. Additionally, the extent 
of overlap also depends on system size. Although the 
longest wavelength modes tend to have plane wave char¬ 
acter, increasing the system size simultaneously generates 
more highly localized modes below the Boson peak, and 
these are most efficient in finding soft spots. These re¬ 
sults call for a systematic study of finite size effects on 
the low energy vibrational spectrum of disordered solids 
and their relationship to rearrangements. 
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